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The dependence of macroscopic detonation properties of a two-dimensional diatomic (AB) molecular system 
on the fundamental properties of the molecule were investigated. This includes examining the detonation veloc- 
ity, reaction zone thickness, and critical width as a function of the exothermicity (Q) of the gas-phase reaction 
(AB — > (1/2)(A2 + B2)) and the gas-phase dissociation energy (D* B ) for AB -> A + B. Following previous 
work, molecular dynamics (MD) simulations with a reactive empirical bond-order (REBO) potential were used 
to characterize the shock-induced response of a diatomic AB molecular solid, which exothermically reacts to 
produce A2 and B2 gaseous products. Non-equilibrium MD simulations reveal that there is a linear dependence 
between the square of the detonation velocity and each of these molecular parameters. The detonation velocities 
were shown to be consistent with the Chapman- Jouguet (CJ) model, demonstrating that these dependencies arise 
from how the Equation of State (EOS) of the products and reactants are affected. Equilibrium MD simulations 
on microcanonical ensembles were used to determine the CJ states for varying Qs and radial distribution func- 
tions to characterize the atomic structure. The character of this material near the CJ conditions was found to be 
rather unusual consisting of polyatomic clusters rather than discrete molecular species. It was also found that 
there was a minimum value of Q and a maximum value of D^ B for which a pseudo-one dimensional detonation 
could not be sustained. The reaction zone of this material was characterized under both equilibrium (CJ) and 
transient (underdriven) conditions. The basic structure is consistent with the Zeldovich- von Neumann-Doring 
model with a sharp shock rise and a reaction zone that extends to 200-300 A. The underdriven systems show a 
build-up process which requires an extensive time to approach equilibrium conditions. The rate stick failure di- 
ameter (critical width in 2D) was also found to depend on Q and B . The dependence on Q could be explained 
in terms of the reaction zone properties. 



I. INTRODUCTION 

The simplest theory of detonation is that of Chapman and 
Jouguet (CJ) Jjl HI Eh- In this one dimensional theory the 
shock rise and reaction are treated as instantaneous. On a 
pressure-specific volume (P-v) state diagram the point of tan- 
gency between a Rayleigh line (an expression of the conser- 
vations of mass and momentum across the detonation front 
traveling at a given velocity) and a Hugoniot (conservation of 
energy) is the CJ state. The slope of the Rayleigh line is pro- 
portional to the negative of the square of the product of the 
initial density (p ) of the material and the detonation velocity 
(u s ). If u s is any slower than the CJ value (it a3 -)> the Rayleigh 
line does not intersect the Hugoniot, and there is no solution to 
the conservation equations. In this light the CJ state is deter- 
mined by the Equation of State (EOS) of the products and the 
initial state of the reactants, and that determines the minimum 
detonation velocity (= u s j) for the system. This hypothesis 
predicts the detonation properties of high performance high 
explosives (HE) reasonably well despite its crude assumption 

A more detailed model is the classical theory of detona- 
tion due to Zeldovich 01, von Neumann [5], and Doling |6] 
(ZND) which, following the initial shock compression, allows 
the molecules of a high explosive (HE) to react and expand. 
This is represented by a pressure profile, whose principal fea- 



tures are ( 1 ) an instantaneous shock rise to a state where the 
reactants are heated and compressed, and this is typically re- 
ferred to as the von Neumann (VN) spike; (2) a fixed-width 
reaction zone, in which reactions provide the chemical energy 
to maintain the detonation wave as density and pressure de- 
crease; and (3) a Taylor wave of the rarefying product gases. 
In the case where the detonation is supported by a driving pis- 
ton, there will be a constant state in the pressure profile from 
some point behind the reaction zone back to the piston. If 
the piston is driven at the particle velocity of the unsupported 
final state (matching the CJ state at the end of the reaction 
zone), there will be no Taylor wave, and only the reaction zone 
will be observed. If the piston is driven at a greater velocity 
than this critical value, then the detonation velocity will be in- 
creased. The detonation is now said to be overdriven, and the 
flow in the constant zone is subsonic in the frame of the deto- 
nation front. For the case where the detonation is underdriven 
with respect to the CJ conditions, it should asymptote to the 
minimum detonation velocity determined by the CJ state, and 
only the Taylor expansion will be affected by a disturbance 
behind the final state ||l|] . 

Molecular Dynamics (MD) simulations are well suited to 
test these and other aspects of detonation theory and their as- 
sociated models under controlled microscopic conditions as 
demonstrated by wo rk g oin g back over a decade H H E3 
d E El 13 E3W El El- With MD it is possi- 
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ble to control the inherent material properties, initial mate- 
rial state, and confinement conditions of the simulation. For 
example, using a predissociative Morse potential, Maffre et 
al. performed a preliminary study of "hot spots", which arise 
at heterogeneities. II 211 . Similar work was also done using 
voids and gaps as the heterogeneities fnl [TsTl . Monte Carlo 
techniques have been incorporated with MD to find thermody- 
namic properties and the Hugoniot of a system of hard spheres 
lU3ll and, more recently, a reactive model close to the one 
used here |9]. Tests of the dependence of the critical flyer 
plate velocity needed to initiate detonation on the flyer plate 
thickness 11511 have been studied. Rice et al. have character- 
ized some aspects of the reaction mechanism 1 10] and demon- 
strated correspondence to hydrodynamic theory for a model 
system |9). Attempts have been made to connect the micro- 
and macroscales with MD and hydrodynamic codes 11911 . 

Several of these simulations ABB damnum m 

[nl CH fl^l have been conducted in two dimensions using a 
Reactive Empirical Bond Order (REBO) potential |7|, rep- 
resenting a simple material composed of two atom types (A 
and B). REBO is a modification of Tersoff 's Empirical Bond 
Order (EBO) potential 1 20] . The restriction to two dimen- 
sions allows significantly greater spatial and time scales to be 
accessed for given computational resources, though three di- 
mensional simulations have also been pursued. The process of 
shock-induced reversible chemical reactivity, converting reac- 
tant AB molecules exothermically into A2 and B2 products, is 
represented in the AB model of Brenner et al. 0. With this, 
it has been demonstrated that non-equlilibrium MD (NEMD) 
simulations using REBO potentials produce detonations con- 
sistent with continuum theory and experimental observations 
lUlll . However, most of these prior studies with REBO have 
focused on using a single set of materials characteristics and 
energetic parameters. It has been found that seemingly sub- 
tle variations in the model or parameters can result in rather 
dramatic changes in behavior. Also, because of computa- 
tional requirements, the simulations have been limited in scale 
such that significant features are not always resolved. Utiliz- 
ing large-scale simulations done with the SPaSM parallel MD 
code 12111 . our aim is to extend the work of Haskins et al. 11611 
by investigating the dependencies of the detonation velocity, 
the reaction zone thickness, and the critical width that a HE 
must have in the transverse directions to sustain detonation on 
these atomistic energetics. This will be done with controlled 
variations in the fundamental microscopic energetic quanti- 
ties, namely the exothermicity (Q) of the reaction and AB 
dissociation energy (D AB ), in order to document more thor- 
oughly the relationship between microscopic properties and 
macroscopic behavior. 

This paper is organized as follows. In Sec.|II|we lay out the 
details of the potential and the simulations used. In Sec. Mil 
using equilibrium microcanonical (NVE) ensembles, we map 
out the product Hugoniot for set values of parameters and 
compare the expected CJ velocities of the detonation fronts 
to those found by unsupported NEMD simulations. The CJ 
states of these materials are also characterized here. 

In Sec. IIVI we study the relationship of Q to the equation 
of state (EOS) in order to better understand the results from 



Sec. [HI] In Sec. [V] we characterize the width of the reaction 
zone and compare this to the critical minimum width (W c ) 
which a 2D sample must have in the direction normal to det- 
onation front's propagation in order for that detonation to be 
sustained. (In 3D cylindrical samples, of course, the critical 
width analog is known as the failure diameter, and the experi- 
mental setup is called a rate stick.) 



H. METHODS 

There are several versions of REBO used in the related lit- 
erature. The version used here is due to Brenner et al. \T\. In 
it the binding energy of an N-atom system takes the form, 

N N 

Eb = ^Y,{fc(r ij )\y R (r ij )-B ij VA(r ij )] + VW^)}, 

i j>i 

(1) 

where ry is the distance from atom i to atom j. Va and Vr 
are the attractive and repulsive terms, respectively, of a Morse 
intramolecular potential Va — Vr, and By = (By + BjCj/2, 
which contains the effective valence and is designed to create 
dimers. Depending on the local environment, By varies from 
to 1. If atom i has no neighbors (defined by f c (r)) other 
than j, By — 1 and the full Morse attraction is felt. On the 
other hand, if i has two neighbors, j and k, and ry < r^fc, 
< Bik < By < 1, i.e., the ij and ik attractions are both 
reduced, but more so for the pair farther apart (ik) than for the 
nearer pair (ij). The effect of this is to introduce a preferential 
valence of one for each atom. A weak intermolecular van 
der Waals (Lennard- Jones form) interaction Kdw stabilizes a 
crystalline AB molecular solid, at least at low temperatures. 
f c is a cutoff function which smoothly takes the potential to 
zero within a finite distance. The rest of the precise functional 
forms and parameters are contained in the errata of |7]. 

Isolated XY molecules (X,Y 6 {A,B}) have binding en- 
ergies Z?g Y (since By = 1); these are the fundamental pa- 
rameters which we will vary from their baseline values |7], 
D AA = D BB = 5.0 eV and D AB = 2.0 eV. These two param- 
eters (constraining D AA = D BB ) are related to the exother- 
micity (Q) through: 



2AB -> A 2 + B 2 + 2Q 



= Dt 



Di 



(2a) 
(2b) 



D AB is just the energy required to dissociate an isolated AB 
molecule: 



AB + D 



AB 



B. 



(3) 



Q is varied from 1.5 eV to 10.0 eV by holding D AB constant 
at 2.0 eV while D AA is varied from 3.5 eV to 12.0 eV. is 
varied from 0.5 to 3.5 eV, by varying D AB and D AA together 
so that their difference, Q, is constant at 3.0 eV. 

In each NEMD setup a stable A 2 flyer plate, four cells thick 
(two A 2 dimers per cell), impacts a 2D metastable AB herring- 
bone lattice at z = and a velocity of 9.8227 km/s (Fig.^. 
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FIG. 1 : (a) A snapshot of a magnified section of the initial sample in 
a flyer driven unsupported NEMD simulation of detonation at t = 0. 
The A atoms are black and the B gray. The A2 flyer plate can be seen 
to the left of the x axis. The box encloses a herringbone lattice cell, 
consisting of 2 AB dimers. 

(b) A snapshot of the detonating sample with the detonation front 
moving toward the right. 



The sample has an initial temperature of 1 1 .6 K. (The poten- 
tial parameters roughly correspond to N2, with a correspond- 
ingly low melting and boiling point.) The resulting unsup- 
ported shock travels to the right, z > 0. In £ the boundaries 
are culled, i.e., particles are lost when they exit the simulation 
cell. In the lateral x direction the boundaries are periodic in 
Sec.|Ill]to study planar detonations or culled, but padded with 
free space, in Sec.[V]to determine failure widths. In compu- 
tations involving the reaction zone thickness and detonation 
velocity, the sample initially is at least 48 lattice cells in x and 
600 in z. The simulation is allowed to run for at least 30.54 ps 
with a timestep dt = 0.25 fs. The minimum duration is the 
same for trials determining W c versus either Q or D^ B . The 
maximum length in cells for these calculations is 360 cells in 
z. In the equilibrium MD calculations presented in Sec. IIVI 
the samples are 25 x 25 cells 2 with periodic boundary condi- 
tions. The simulations are run for 40 ps with measurements 
averaged over the last 30 ps. 



III. DETONATION VELOCITY AND THE CJ 
CONDITIONS 

Haskins et al. report a linear dependence of u 2 on Q for a 
similar AB material 1 16]. We repeat this study with a slightly 
different REBO potential and further it by varying the AB dis- 
sociation energy (D^ B ). One expects the velocity to increase 
with Q because the increased exothermicity of the reaction in- 
creases the temperature and pressure of the products. Fig. [5] 
confirms the linear relationship between u 2 and Q for this par- 
ticular system. The differences between values in Fig.|2]and 
lUrjll are due to differences in other REBO parameters and the 
flyer's thickness and impact velocity. Here, impact velocity 
is not raised above 9.8227 km/s in order to avoid the "fast 
detonation regime" 111 211 (which may have been an artifact of 




FIG. 2: Square of the detonation velocity (tt s ) vs. AB dissociation 
energy (D^ B ) and exothermicity (Q). The flyer thickness and im- 
pact velocity are held fixed, so detonation is not sustained in regions 
where it is reported to in 1 1 611 . Failure is indicated by the shaded 
boxes. A sustained detonation is defined as not failing before the 
right edge of the sample. Error bars are smaller than the size of the 
symbols. The scale on the right gives the kinetic energy of an atom 
traveling at u s . The filled symbols are from NEMD simulations with 
a free boundary at z — 0. The empty symbols are from NEMD 
simulations with a momentum mirror at z = 0. The linear fits are 

57.251 + 16.694D* 8 . 



u 2 s = 33.205 + 19.137Q and ul 



the model used there). At values of Q < 1.5 eV, the lin- 
ear relationship begins to fail, and the system would not sus- 
tain a propagating detonation for Q = 1.3 eV. This failure 
point likely arises because the reaction rate (determined by 
the temperature at the initial shock front) has become suffi- 
ciently slow such that it does not approach completion within 
the subsonic region of the reaction zone 1 1]. The temperature 
at the shock front can be estimated from the kinetic energy of 
the shock front, given by the right-hand axis labels in Figure|2] 

The dependence of u 2 on D^ B is also found to be linear and 
increasing (Fig.|2jl. Our initial expectation was that the vari- 
ations in D^ B would primarily affect the activation energy of 
the reaction and could quench the detonation when the acti- 
vation energy became too great to be readily overcome at the 
temperature of the initial shock state. The failure to main- 
tain a propagating detonation for D^ B > 3.7 eV is possibly 
a manifestation of this. The strong dependence of u 2 on D^ B 
indicates that other aspects of the system are likely being af- 
fected by this perturbation. 

To understand these relationships more thoroughly, we turn 
first to the basic test of standard detonation theory, which is 
to compare predictions based on the CJ state determined from 
equilibrium MD simulations with the evaluation of detonation 
propagation from NEMD studies. The theoretical CJ state for 
a ID detonation is at a sonic point at thermoc hemic al equilib- 
rium |l|]. If this is an improper assumption and the reaction 
is incomplete at the arrival time of the sonic point, it could 
account for a discrepancy between the theoretical and actual 
detonation velocities. For a slightly different REBO potential, 
Rice et al. found a detonation velocity from an unsupported 
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simulation that was 6.1% lower than that found from the cal- 
culation of their CJ state |9]. 

Our procedure for locating the CJ state is described here 
and similar to that of Rice et al. |9J]. At different values of 
specific volume (v), sets of microcanonical (NVE) MD simu- 
lations of 2500 particles are run for 40 ps, providing sufficient 
time to equilibrate. The onset of equilibrium is determined 
by the shape of the time evolution of the average properties 
of simulation. After each of these reaches a plateau, as de- 
termined by visual inspection, the simulation is allowed to 
continue. A runs-above-the-median test is performed to de- 
termine that the curves are fiat with only random fluctuations. 
At each value of v, the value of specific internal energy (E) is 
sought that is a solution of the Hugoniot jump condition, 

^(v -v)P = E-E (4) 

(energy conservation in which Po = 0), where P = \ (P zz + 
P xx ) is the hydrostatic pressure, P a p is the negative of the 
corresponding component of the stress tensor and has ideal 
and virial components, and the subscript indicates the state 
in front of the detonation front. Once (E) is determined for 
the present value of v, NVE ensemble averages (x) of other 
thermodynamic quantities are computed by linear interpola- 
tion. 

By repeating this procedure for several values of v, the 
product P-v Hugoniot is determined (Fig. Using the re- 
maining jump conditions, mass conservation: 

u p /u s = (v Q - v)/v (5) 

and momentum conservation in which Po = 0, 

u s u p = v P, (6) 

one can find u s as a function of the particle velocity at the final 
state (u p ). The minimum possible value of u s is the CJ value 
(jl, so the minimum of u s versus any thermodynamic parame- 
ter is at the CJ value of that parameter (e.g., see Fig.|4}. We use 
the minimum determined by these means as iterative approxi- 
mations of the CJ state (see Table|I]under the CJ-Interpolation 
column). We refine the process by fitting a quadratic through 
the points surrounding the minimum. We must be careful to 
include a domain small enough that a quadratic is a good ap- 
proximation to the data, yet large enough to include more than 
three data points in order to get a proper estimate of the error 
from the goodness-of-fit parameter. 

To test the CJ results found from these NVE simulations, 
we run a supported detonation simulation in which an in- 
finitely massive driving piston impacts the AB sample at the 
CJ-determined velocity (u p j). This should establish a con- 
stant zone behind the front that should be at the CJ state, 
with the front propagating at the CJ conditions. Measure- 
ments are taken from this run and can be found in Table U 
under the Supported-NEMD column. From Table U one can 
see that many of the values fall within error of one another, 
although there are some slight discrepancies. It should be 
noted that the supported-NEMD simulations can include tran- 
sients from the initiation and/or build-up processes |2j], and 
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FIG. 3: A P-v state diagram of the equilibrium product Hugoniot. 
The solid curve is a guide to the eye. The dotted/dashed curves are 
Rayleigh lines plotted using the initial state (Po = 0.0, v/vo = 1.0) 
and a slope of -ulpo, where p is the mass density and u s for each 
curve, from the steepest down, is the median value of shock velocity 
taken from the supported detonation, the EOS calculation, and the 
unsupported detonation. The arrows represent a series of simulations 
at constant v used to find the datum to which they point. The box is 
a magnification. 



these will be explored further below. Still, this discrep- 
ancy of 0.35% difference from the unsupported detonation of 
u s = 9.70961 ± 0.00054 km/s is quite small, especially since 
previous evaluations of this same system gave larger discrep- 
ancies (9.3 km/s in |7|] and 9.5 km/s in |8]). This is also less 
of an error than Rice et al. found for their variation. 

We repeated several of the tasks performed on Q = 3.0 eV 
for several more values of Q. In Fig. [5] the equilibrium P- 
v Hugoniots for several values of Q are shown. Notice for 
Q = 6.0 eV that the curve has a nice hyperbolic shape, for 
which the CJ state for this Q is easy to determine. From it, as 
was done above, we find u S j to be 12.20552 ± 0.0020 km/s. 
The unsupported NEMD simulation with the free boundary 
gives 12.2031 ± 0.0017km/s, a difference that is smaller than 
the error bars. As Q decreases, the equilibrium Hugoniots 
flatten out, and it becomes more difficult to identify the exact 
point of tangency. At Q = 2.0 eV the Hugoniot is well rep- 
resented by a linear fit. The Hugoniot for Q = 1.5 eV has a 
negative curvature (convex) section indicating the presence of 
some sort of phase transition (Fig.|5Jl. It should be noted that 
Q = 1.5 eV is close to Q c , the value for which a detonation 
can not be sustained. 

A comparison of the predicted and measured values of u s 
for different values of Q are given Table [H] It is found that 
there is very good agreement between these values, with no 
discernible discrepancy for Q = 6.0 eV, and only a moder- 
ate discrepancy of 2% for Q = 1.5 eV which is very near the 
failure point. The latter point supports the idea that failure 
is occurring because of inadequate completion of the reaction 
when the sonic point passes. Overall though, it can be con- 
cluded that these systems are responding in accord with the 
simple predictions of the Chapman- Jouguet hypothesis. The 
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TABLE I: Values measured or defined in a process which determines the CJ state through a series of microcanonical equilibrium simulations 
or in the constant zone of a non-equlilibrium simulation of a critically supported detonation (Sec. HVi . v is the specific volume. u s is the shock 
velocity. u v is the particle velocity at the final state. T is the absolute temperature. U is the potential energy. E is the specific internal energy, 
and P is the 2D pressure. Values marked with asterisks are input into the corresponding simulation(s). The parentheses indicate the error in 
the last two digits of the corresponding reported value. 





CJ Interpolation 


Supported NEMD 


Q(eV) 


1.5* 


3.0* 


6.0* 


3.0* 


v/vo 


0.50182(24) 


0.57164(81) 


0.6171(19) 


0.5747(16) 


u s (km/s) 


8.0720(11) 


9.6758(90) 


12.2055(20) 


9.7360(40) 


Up (km/s) 


4.0212(13) 


4.1446(40) 


4.673(22) 


4.144* 


(k B T) (eV) 


0.4017(14) 


0.7839(43) 


1.7550(90) 


0.7791(56) 


(U) (eV) 


-0.25156(79) 


-0.5605(21) 


-1.19657(90) 


-0.5746(72) 


(E) (eV) 


0.15012(65) 


0.2234(22) 


0.5584(84) 


0.2234306(15) 


P (eV/A 2 ) 


0.70591(33) 


0.8726(13) 


1.2411(62) 


0.8541(83) 





FIG. 4: Shock velocity (u s ) vs. particle velocity (u p ). Similar plots 
can be made for other thermodynamic variables on the abscissa. Each 
point corresponds to a different value of specific volume (v), which 
increases along the abscissa in the direction opposite that of u p . A 
constant line is drawn at the determined minimum value of u 3 . The 
dotted part of the curve indicates that which is not dynamically ac- 
cessible. Instead NEMD simulations with pistons moving at u p less 
than the value at which the minimum is located ideally should follow 
the constant line. The open circles show the results of piston driven 
NEMD simulations where the piston is driven at u p . 



variations in u s with the molecular parameters is therefore de- 
termined by how these effect the Equation of State (EOS) of 
the products. These aspects will now be examined in further 
detail. 

The curves in Figs. [5] and [4] are rather unusual. For a mate- 
rial behaving like an ideal gas, the former would be expected 
to have a hyperbolic shape and for the latter a parabolic shape. 
Here, for large Q (= 6 eV) the resulting curves do have that 
type of behavior. However, as Q is reduced, the Hugoniot flat- 
tens until it eventually gains a convex section, and the u s -u p 
plot evolves into a curve with a double minimum. Therefore 
we examine the CJ states more closely by simulating a mi- 
crocanonical ensemble at the determined CJ values of density 



FIG. 5: Equilibrium Hugoniots for several values of Q. The CJ state 
is determined for the cases in which Q = 6.0, 3.0, and 1.5 eV. The 
curve for Q = 2.0 eV is well fit by a straight line of the form 
P = 1.69(1 — v/vq). The dotted lines are all guides to the eye. 
As Q decreases, the Hugoniots go from positive curvature to a zero 
curvature for Q — 2.0. The curve for Q = 1.5 eV contains a section 
of negative curvature, indicating a phase transition. The arrows point 
to the determined CJ states. 



TABLE II: As Q is increased, the predicted (u S j) and measured val- 
ues (u s ) of the detonation velocity agree more. It is hypothesized 
that the increased exothermicity causes the reaction to come closer 
to completion by the sonic point. 



Q(eV) 


u S j (km/s) 


m 3 unsupported 


% difference 


1.5 


8.0720(11) 


7.913(16) 


« 2.0 


3.0 


9.6758(90) 


9.70961(54) 


« 0.35 


6.0 


12.2055(20) 


12.2032(17) 


« 0.0 



and internal energy. 

A snapshot of the NVE-at-CJ simulation for Q = 3.0 eV 
is shown in Fig. [6] along with the corresponding radial dis- 
tribution functions (RDF), g(r), for particles of the same and 
different type. It is evident from Fig. [6] that, at the CJ state, 
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the system is dominated by short AA and BB contacts, which 
would be expected for the product species. The interatomic 
distances for these (w 1.1 r e ) is slightly larger than that de- 
fined for the isolated molecules (1.0 r e ). By the form of the 
potential, the introduction of a third particle within the de- 
fined bonding distance weakens the attractive part of the po- 
tential and moves the minimum of the bonding potential out 
to greater distances. This aspect is readily explained by the 
high density at CJ. 

What is more intriguing is that there is also a significant 
number of close A-B interactions at w 1.2 r e , showing that 
the system has not evolved into a simple mixture of A2 and 
B2. Also, there is a second peak in both RDFs at around 
r/r e = 2.0. This second peak indicates that there are clusters 
of atoms forming at this compression. Looking at the snap- 
shot embedded in Fig. [6] one can see that there appears to be 
linear oligomers forming. 

For comparison Fig.0shows the RDF of an NVE simula- 
tion of AB at a temperature above its melting point and at the 
initial density, v = Vq. It has a strong peak for the AB dimers 
at r/r e — 1 as expected. Beyond this, it was expected that 
there would be a broad peak at 2.8 r e , which would be the 
van der Waals minimum. Although this was observed, some- 
what surprisingly, there is a sharp peak just above r/r e = 2 
which is present for all three atomic combinations, as well as 
another peak at 3 r e . The former peak is caused by the posi- 
tive slope of a section of the inner cutoff spline for the V v dw 
term, which creates an additional minimum in the multibody 
potential. Examination of the inset show that the molecules 
still tend to line up, although the molecules are clearly sep- 
arate diatomic species. The peak at r/r e — 3 is probably 
a result of this alignment and is the distance from a bonded 
partner of one atom to the non-bonded neighbor which is at 2 
r e . These features highlight some of the unusual aspects that 
can arise from these complex interaction potentials. 

The notable difference between the RDFs is that the one 
for the melted AB goes to zero between r/r e — 1.2 and 1.8, 
which highlights the clear diatomic nature of that system. For 
the CJ state simulation, there is substantial intensity all across 
that region, and the RDF value barely drops below 1 . Since 
there are no dissociated atoms or oligomers in the melted AB, 
it is reasonable that we see a domain above r/r e = 1 in which 
a particle will not have a neighbor. Any third particle and its 
bound partner will be repelled by virtue of the Va term if they 
approach much closer than 2 r e to a particle in another dimer. 
At the CJ state density, this repulsion breaks down and clusters 
of particles form. Similar behavior has also been observed in 
the systems that Rice et al. studied, where describing the CJ 
state as one of atomization is a fair description. l9l ll0ll . 

To examine this behavior more closely, we examine the 
RDFs along the Hugoniot. In particular, we consider the case 
of Q — 1.5 eV which has a convex section that indicates a 
phase transition. The RDFs for several states along that Hugo- 
niot including v/vo = 0.64, the volume at which the phase 
transition occurs, are shown in the Fig.[H)| For the CJ state, 
v/v n = .502, the RDF is similar to that for the CJ state at 
Q = 3.0 eV. The first peak occurs at 1.2 r e with no strong 
minimums occurring at longer distances. At v/vq = 0.64, 



TABLE III: The minimum peak value of the barrier to be overcome 
when converting AB to BB or AA in 1 dimension (-B a ,o, where the 
second subscript indicates the angle between AB and BB) for several 
values of D^ B . Detonation cannot be initiated at D^ B = 4.0 eV for 
the given flyer thickness and velocity. All of the measurements of 
E a fl have an error of ±0.0125 eV. 



(eV) 


E a ,0 


3.0 


0.1875 


3.25 


0.2375 


3.5 


0.2625 


3.75 


0.2875 


4.0 


0.3375 



a deep minimum develops at sa 1.5 r e along with peaks at 
r/r e = 2.0 and 3.0. From the snapshot at v/vq — 0.64 (See 
inset.), one may be convinced that both distinct dimers and 
extended oligomers are present. At v/v = 0.75, the peak 
sharpening and minimum development are more distinct, and 
the profile is more reminiscent of the AB melt illustrated in 
Fig.0 

To investigate the cause of these clusters, we plot potential 
energy surfaces (PES) for the interatomic distances of three 
inline particles as Rice et al. | 10] did for a different version 
of REBO. From Fig. [8] which represents an A-A-A configu- 
ration, one can see that, at a distance of about 1.35 A, there 
is a minimum that allows for trimer formation. This config- 
uration does not exhaust the possible configurations that may 
occur during reaction because it is constrained to be a linear 
conformation. The interaction with neighboring atoms, par- 
ticularly at CJ-type conditions, is also ignored and this could 
alter the absolute value of the activation energy (E a ) signifi- 
cantly. However, as Tablelllll suggests. E a seems to increase 
mono tonic ally with D^ B . This would rationalize that as D^ B 
is increased the material eventually fails to detonate given the 
same initiator strength. 



IV. EXOTHERMICITY'S RELATIONSHIP TO EOS 

We now turn to understanding the linear relationship be- 
tween u 2 and Q. One simple model is suggested by Fickett 
et al. OIL wno derive a linear relationship between Q and u 2 
by adding Q to the incomplete equation of state of the initial 
state of a polytropic gas, an ideal gas with a constant specific 
heat. An expression for the change in specific internal energy 
becomes E — Eq = (Pv — PqVq)/(^ — 1) — Q, where 7 is 
the adiabatic gamma |22]. One can then eliminate E — Ea 
with the Hugoniot jump condition (Eq. [4}, thus solving for 
P. It is assumed that Pq = 0. One can use the condition 
that Rayleigh line is tangent to the Hugoniot at the CJ state to 
arrive at u 2 = 2(^ 2 — l)Q/m, where m is the mass of the 
reactants. If it is assumed that this EOS accurately describes 
our potential, 7 « 2. A typical conventional HE has a 7 « 3 
01. Since the volume at the CJ state is given by the expres- 
sion v/vq — 7/(7+ 1), this relationship does rationalize the 
somewhat greater compression observed here (v/vq < 0.67) 
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FIG. 6: Radial distribution function for the CJ state for Q = 3.0 eV. 
(?ab measures the chance of finding particles of opposite types a dis- 
tance r apart divided by the probability if the atoms were randomly 
distributed. #aa/bb indicates the probability of finding particles of the 
same type a distance r apart, again, divided by the probability if the 
distribution were random. The inset shows a section of a snapshot of 
the simulation. A atoms are black and B are gray. 




FIG. 7: RDF for melted AB. The first peak, the maximum of which is 
about 14, indicates the presence of dimers. The peak above r jr e = 2 
is due to the inner cutoff spline in the V v aw term of Eq.Q The peak 
at r/r e = 3 is probably caused by the arrangement of dimers. 
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r *-4 <A) = 



FIG. 8: Potential energy in eV vs. separation distances in A of three 
inline A (or B) atoms. The "*" marks the contour labeled by the 
nearby number. Notice the local minimum at (1.3, 1.3). It allows for 
trimer formation. No positive contours are shows, as would fill in 
the section near (1,1). The case for inline A-B-A (or B-A-B) looks 
similar yet shallower. 




FIG. 9: Potential energy vs. separation distances of an A-A-B (or 
B-B-A) configuration. 



compared to that for conventional HE's (v/vq ~ 0.75). (It 
should be noted that the values of v/vo given in Table H] are 
smaller yet and imply a value of 7 « 1.) However, as we have 
demonstrated that the product's EOS does not behave like an 
ideal gas, this does not serve as a good basis to explain the 
observed relationship between u 2 s and Q. 

An alternative explanation is suggested by an examination 
of the Hugoniot curves shown in Fig.|5] There, in the region of 
0.55 < v/vq < 0.7 (which spans the region containing the CJ 
states), there appears to be an approximately linear offset of 
the curves. This suggests that a Mie-Gruneisen EOS with an 
unspecified reference Hugoniot might be suitable, where it is 



assumed that the value of T/v is dependent on v. Truncating 
the Taylor expansion around a reference Hugoniot as in 1 22], 
we get 



P « P R + T(E - E R )/v, 
where the Griineisen gamma 

dP 
dE 



(7) 



(8) 



We substitute Eq.|4]into Eq.Qand, because we take E as rela- 
tive to the products, set Eq = q, where q is the specific heat of 
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vlv = 0.502 (CJ) 

= 0.64 
vlv„ = 0.75 




FIG. 10: Same-type radial distribution function for an exothermicity 
Q — 1.5 eV at states along the Hugoniot at volumes shown. Notice 
the maximum at radial distance r/r e ~ 2.75 for the CJ state. It rep- 
resents the van der Waals equilibrium distance. The local minimum 
of the other curves at that point are still greater than unity. Notice 
the local maxima for the other two curves at r/r e — 3.0. They sug- 
gest tetramer formation or dimer alignment as with the AB melt. The 
inset is a snapshot of the simulation for v/vo = 0.64. 



reaction. Upon rearrangement we get for the P-v Hugoniot 
1 



H 



V (l 



r v 



[vP R + Y(q - E R )] . (9) 



2 v J 

Rearranging Eqs.|5]and|6j we find 



2 D V 

u s = P H - 



1 



This yields the general expression 



(10) 



u 2 s = A(v,T)P R + B(v,T)(q-E R ). (11) 

It should be noted that q does not necessarily equal Q/m. In 
the case where the bond order parameter Bij = 1.0 for the 
initial state, Q /m should be q less the small contribution of 
the van der Waals interaction. When we run an NVE simu- 
lation of 1250 AB molecules at the initial state used in the 
NEMD simulations, we find a total internal energy of -2548.5 
eV. Where Davis measures the zero of energy as cold products 
ll22ll . our calculations use a zero of cold dissociated atoms. 
Our measurement here finds Eq « D^ B and is consistent with 
Bij = 1.0 and a 2% van der Waals contribution of a few 
neighbors, therefore Q/m w q to very good order, and the 
two can be used interchangeably. 

The validity of this Mie-Griineisen approximation can be 
tested by the linear dependence of u 2 on Q for constant vol- 
umes. A few of these plots are given in Fig. ^3 which shows 
that good linear relationships are found for the three selected 
volumes. This shows that the Mie-Griineisen EOS is a good 
approximation for our model for low v (high compressions). 
Since the data were generated during the process of seeking 
the minimum u s for each Q, interpolation was used to find 



values at common abscissas. The coefficients of the fits gen- 
erated by the method in Fig.[^ are then plotted in Fig.ll2lvs. 
v/vq. Here, A' is the y-intercept of those linear fits, and B is 
the slope. The value A' is distinct from A in Ea.llllbecause 
it now includes a contribution from E R . However, the value 
of B is exactly the same as defined in that equation. The fact 
that all of the lines in Fig.[^cross at Q = 2.0 eV leads one 
to conclude that A' and B are linearly dependent. The third 
curve on Fig. ^]is an attempt to show this for the relation 
A' = 77.5 — 2.0S. 

The net result of this analysis is that, over the range of 
0.55 < v/vq < 0.7, B is reasonably constant with a value 
of 17-19 (km/s) 2 / eV. This result, inserted back into Eq.ll II 
is then in good agreement with the initial observation shown 
in Fig. |2j that the slope of the dependence of u 2 on Q is 
pa 19 (km/s) 2 /eV. The ability to approximate the product 
EOSs as linear Mie-Griineisen offsets of one another, at least 
in the region near the CJ states, clarifies the origin of this lin- 
ear dependence for this system. 

Using Eqs.[9]E3 andll II we can then solve for T to find 



2B 



1 _ JL JL 

Vo J Vq 



2 + B 1- 



J» 



and plot the result in Fig.^] Solving for B we have 

r 



(12) 



B(v) 



1 - 



1 - 



M)) 



(13) 



Note that (1 — v/vq) is the compression (77). From these equa- 
tions, it can be seen that a constant value of B does not specify 
a constant value of V, but rather a particular dependence of T 
on v. Conversely, a constant value of V would also specify a 
particular dependence of B on v. However, from these two 
figures it is observed that both B and T are rather weak func- 
tions of v in the CJ region. This aspect of the EOS behavior 
(which we are unable to directly adjust) is the origin of the 
linear dependences observed in Fig. [2] 



V. REACTION ZONE THICKNESS AND RATE STICK 
FAILURE DIAMETER 

Having established the steady-state properties of these sys- 
tems, we now turn to the transient properties related to the 
reaction zone thickness. There are two means that could be 
employed to determine the character of the reaction zone for 
these simulations. The most direct means is to drive the sys- 
tem with a piston whose velocity is matched to that of the CJ 
conditions. In this case, the one dimensional profile should 
exhibit a rapid shock rise up to the von Neumann state, fol- 
lowed by a relaxation to the CJ state, which is then a constant 
zone that extends back to the piston. There will be some tran- 
sients present because the initiation process incurs a slight de- 
lay so that an equilibrium state is not immediately established, 
but otherwise this is a direct method if the CJ conditions are 
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FIG. 11: Square of the shock velocity (ttg) vs. exothermicity (Q) for 
select constant values of specific volume (v). The lines are linear 
fits. Notice that the fit is better for the higher values of v. Notice, 
also, that the lines cross at v/vq = 2.0 eV. This indicates a linear 
dependence between the slopes and the y-intercepts. 




FIG. 12: Coefficients of the fits as described in Fig. llll vs. specific 
volume for D^ B =2.0 eV. A' is the y-intercept and, as in Eg. II II B 
is the slope. The dashed line is the calculation shown in the legend. 



known. Less direct is to study a system with either an un- 
derdriven piston or perhaps initiated with a flyer plate. Such 
systems will initially propagate at less than the CJ conditions 
because the release waves can eat into the reaction zone. As 
the detonation proceeds and the release Taylor wave becomes 
more spread out (approximating a more steady condition be- 
hind the reaction zone), these will eventually "build-up" to 
a steady detonation. It can be difficult to determine exactly 
when the system has evolved to a steady-state condition. 

A comparison of these two approaches is shown in Fig. 1171 
For the case of a piston matched to the CJ conditions, a steady 
solution evolved quite rapidly as determined by both the det- 
onation velocity and the constant properties of the following 
zone. The reaction zone, defined to be the distance from the 
shock front to the point where the transient properties are not 
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FIG. 13: Griineisen T vs. specific volume for D^ B — 2.0 eV. 



distinguishable from the thermal fluctuations of the constant 
zone, extends out to f» 300 A. A somewhat better character- 
ization is probably the distance at which the particle veloc- 
ity has decreased to 1/e of the initial overshoot. This occurs 
at « 60 A behind the shock front, or a characteristic time 
constant of 0.7 ps. For the transient simulations, it is appar- 
ent that the system is continuing to evolve even after prop- 
agating for 100 ps. The reaction zone length (emphasized 
in the inset) is clearly extending past 100 A. (The unusual 
structure after the reaction zone is likely due to the unusual 
chemical structure of the CJ state, and could be character- 
ized as a phase transition back to the diatomic A2 and B2 
products.) These measurements for the build up time, reac- 
tion zone width, and the detonation velocity are larger than 
those proposed by White et al. |8] for the same model. They 
include in their calculations a time domain in which the deto- 
nation is still subtly building. We, in fact, measure a detona- 
tion velocity (9.5580 ± 0.0013 km/s) close to theirs (9.5 km/s) 
for the unsupported case if we include the front positions at 
10 ps < t < 20 ps in the linear least-squares fitting that de- 
termines u s . This highlights the difficulties in using this ap- 
proach of unsupported detonations. 

It should be noted that the transient system exhibits the clas- 
sic build-up behavior quite nicely. In addition to the extension 
of the reaction zone out to its equilibrium value, there is also 
the gradual increase in the initial shock state up to its equi- 
librium value. As this sets the initial temperature and rate for 
the reaction chemistry, it emphasizes the significance of these 
transient phenomena. The shock velocity is also a sensitive in- 
dicator of the build-up process. While the deprecation of the 
velocity below the steady-state value is not great nor readily 
obvious, it can be discerned by careful observation. Previ- 
ously reported discrepancies between expected and observed 
CJ performance may well have been due to systems not hav- 
ing fully attained steady-state conditions. 

In 12311 Davis shows that surface effects adversely affect 
detonation's ability to propagate because rarefaction waves 
from the side release erode the reaction zone and its support 
structure. For cylindrical charges, this is characterized as a 
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FIG. 14: Rate stick at super-critical width. 



FIG. 15: Three sequential snapshots of a rate stick at sub-critical 
width. Here the particles are colored by bond type. Gray is unre- 
acted, and black is reacted. Rarefaction waves from the rate stick's 
edge erode at the reaction zone thus quenching the detonation. 



failure diameter below which detonation cannot propagate. 
White et al. report that, when the periodic boundaries are re- 
moved from the sample so that the particles are free to escape 
the system, there is a critical width (W c ) that the sample must 
exceed in x, the direction perpendicular to the propagation of 
the shock front, in order that the detonation be sustained. If 
the sample is too thin, rarefaction waves quench the detona- 
tion n. 

W c 's dependence on Q llal and B can be seen in Fig. 1161 
Unlike in II 611. a curve of the form y = a/(x — Q c ), not ex- 
ponential, is fit to the data for Q. Simpler equations fit better, 
but this form is inspired by Hubbard and Johnson's l24ll de- 




fi.O (eV) 



FIG. 16: Critical width (W c ) vs. AB dissociation energy (D ^ ) and 
exothermicity (Q). The upper dash of the error bar indicates the 
thinnest sample to sustain detonation. The lower error bar indicates 
the thickest sample to fail within 360 lattice spaces. The resolution is 
10 A. The lines are guides to the eye with empirical functional forms 
discussed in the text. 



pendence of delay time (td) (described below) on Q and to be 
asymptotic at the point where detonation is known to fail for 
an infinitely thick sample (Q c ). For D^ B a curve of the form 
y = a(e bx — l)/(x — D™?) is fit to the data for similar reasons. 
If Q is raised, the reaction should be more difficult to quench 
because there is more energy available to break neighboring 
AB bonds. The reaction is expected to be faster and the reac- 
tion zone shorter because of the higher temperature resulting 
from the greater energy release. As Q lowers toward Q c , no 
matter how wide the sample is, detonation will not be sus- 
tained. On the other hand, if D^ B is lowered, it is harder for 
rarefaction waves to quench the detonation front because, with 
a lower dissociation energy, a chemical reaction occurs more 
readily. W c is, thus, lowered because detonation is more eas- 
ily sustained. As it rises, D^ B reaches a critical value (D^ B ) 
(probably dependent on Q, flyer thickness and velocity, etc.) 
above which detonation cannot be sustained, no matter how 
wide the sample. 

In a sample with culled boundaries, if b is lowered, it 
is harder for rarefaction waves to quench the detonation front 
because, with a lower dissociation energy, a chemical reaction 
occurs more readily. W c is, thus, lowered because detonation 
is more easily sustained. As it rises, D^ B reaches a critical 
value (D^ B ) (probably dependent on Q, flyer thickness and 
velocity, etc.) above which detonation cannot be sustained, no 
matter how wide the sample. 

Using the idea of delay time (td), during which a reactant 
must remain above a certain temperature in order for a reac- 
tion to occur, Hubbard et al. give an example of how a deto- 
nation's ability to be initiated or sustained may depend on Q 
and D^ B . We use this in lieu of finding reaction rate data. The 
longer this time, the less likely the reaction. In their example, 
td has an inverse dependence on Q and an inverse and expo- 
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nential dependence on activation energy, E a (Eci.ll4> |2 

REl 



(a) 



td 



c v vE a Q 



exp 



REq 



(14) 



where v is the collision rate, R is the molar gas constant, c v 
is the specific heat at constant volume, and Eq is the internal 
energy behind the shock. 

In Rosing and Chariton theory, when 8 = t r , d — W c , 
where is scattering time, t r is the reaction time, and d is the 
diameter of the sample. The rate stick must be thick enough 
that rarefaction waves, traveling at c from the surface, cannot 
sufficiently penetrate the HE in order to scatter the reaction 
zone and quench the reaction in a time equal to the reaction 
time, a = t r (u s — u), where a is the width of the chemical 
reaction zone, u is the average particle velocity in z within the 
zone, and u s is the velocity of the detonation front. = d/2c, 
where c is the speed of sound. Replacing t r with 0, one finds 



W c — d = 2ca/(u s - u), 



(15) 



where u s is a function of d 111 lll23l 12511 . If td oc t r and E a oc 
D^ B , we find that W c oc td, supporting our selection of curve 
fits. 

To check the validity of Eq. [21 we use for a the charac- 
teristic decay length of the reaction zone (60 A) determined 
above. If the ratio 2c/ (u s — u) ~ 2, this would be excellent 
agreement. As the value of u s — u should be approximately 
equal to the local sound speed, this estimation is then quite 
good. Given the approximations involved in making this as- 
sessment, this is certainly a fortuitous agreement, but can be 
taken as support for this analysis. Overall, a more thorough 
understanding of these phenomena is evolving through these 
studies. 



VI. CONCLUSION 

We have analyzed the relationships between the energetic 
chemical properties of a simple, but well studied, model high 
explosive and the physical properties of its detonation front 
through explicit Molecular Dynamics (MD). This provides 
an excellent means to compare to continuum analyses and 
avoids the complication of a multiphase equation of state 
(EOS) and other approximations required in those approaches. 
Most previous MD studies of high explosives have focused 
on the physical properties of the HE as independent param- 
eters (e.g. voids and gaps) to validate this general approach. 
Our analysis has been conducted in order to establish a more 
quantitative comparison to the relevant theories by, e.g., seek- 
ing the CJ state via NVE simulations and showing that it ac- 
curately predicts detonation properties such as the detonation 
velocity (u s ). Despite spurious model anomalies, e.g., a CJ 
state consisting of polyatomic species and unusual proper- 
ties of the product Hugoniot curves, we have observed gen- 
eral consistency between the MD approach and the thermo- 
dynamic equations. Given that the fundamental CJ and ZND 
models pose minimal constraints on the behavior of chemi- 
cal system, this result is not extremely surprising. However, 
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FIG. 17: (a) 129 profiles of the particle velocity on the z direction 
{u z ) separated in time by 50.9 fs are overlapped onto their center 
profile (whose time is indicated in the legend) such that their fronts 
line up. They are averaged and those averages are similarly over- 
lapped. The simulation represented has the default parameterization 
for the REBO potential |7J. It is a flyer driven unsupported simula- 
tion as described in Sec.Ql] Found for this specific parameterization 
in Sec.lTVl u P j is indicated by a constant line segment. The inset is 
a magnification of the region containing the von Neumann spike and 
the CJ state. 

(b) Similar profiles for a critically supported detonation. This simu- 
lation is shorter and thinner than in (a), but it suggests that in (a) still 
has some building to do for the profiles do not seem to settle down to 
Um until about 300 A behind the front. 



we have demonstrated that it is possible to gain a quantita- 
tive understanding of these relationships even for these some- 
what unusual systems. For the analysis of the behavior of u s 
and the critical width as functions of energetic parameters, we 
have explained our results in the context of available theories. 
In the case of u 2 s as a function of the dissociation energy of 
AB, we have not directly determined the reason for its lin- 
ear and increasing relationship, but assume that it also deter- 
mined by the induced EOS changes. In characterizing the re- 
action zone width, we demonstrated that transient simulations 
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must be scrutinized rather carefully in determining the final 
state. We are currently investigating the sources of the model 
anomalies that we have identified, and potential model revi- 
sions that will correct these are expected to result in a system 
that behaves more similarly to conventional high explosives. 
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